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I.  INTRODUCTION 


A  new  theoretical  procedure  is  developed  for  calculating  light  scattering  by  a  spherical 
particle  that  is  on  or  near  a  conducting  plane  surface  (minor).  This  computational  capability  is 
useful  for  both  scientific  investigations  and  engineering  design  studies.  The  scattered  light,  for 
example,  can  be  used  to  detect  and  measure  particles  on  a  flat  reflecting  surface.  This  is  relevant  for 
certain  manufacturing  processes  where  surfaces  must  be  maintained  in  an  uncontaminated 
condition.1  The  present  work  originated  in  the  problem  of  modeling  stray  radiation  in  an  optical 
system  due  to  light  that  is  scattered  from  particles  contaminating  the  surfaces  of  mirrors.  Such  stray 
radiation  can  seriously  degrade  performance  of  the  system.2  In  addition,  if  the  radiation  is  intense, 
the  energy  absorbed  by  the  particle  can  cause  local  heating,  which  may  damage  the  mirror  or 
otherwise  cause  problems.  Another  possible  application  of  the  theory  is  to  the  study  of 
morphology  dependent  resonances  (MDR)  of  small  particles.3  In  recent  years  many  studies  have 
been  carried  out  that  exploit  resonance  effects  in  isolated  particles.4*5  Interesting  new  resonance 
effects  may  be  possible  if  the  particle  is  on  or  near  a  conducting  surface. 

Previous  work2*6*7  on  light  scattering  from  a  particle  on  a  mirror  has  included  born 
experimental  and  theoretical  studies.  Prior  theoretical  calculations  have  been  based  on 
approximation  methods.  The  method  developed  in  this  report  is  not  an  approximation,  but  rather  it 
provides  a  numerically  exact  solution  to  Maxwell's  equations. 

Section  II  begins  with  a  general  discussion  of  the  scattering  problem  and  a  derivation  of 
the  method  of  images.  This  method  allows  one  to  replace  the  system  consisting  of  the  sphere  and 
the  conducting  plane  with  an  equivalent  system  consisting  of  the  sphere  and  an  image  sphere.  The 
equivalent  two-paiticle  system  is  then  solved  by  the  multipole  expansion  method.8*9  Mirror 
symmetry  is  exploited  to  simplify  the  calculation.  The  section  concludes  with  a  derivation  of 
formulas  for  the  extinction  and  absorption  cross  sections.  Except  for  the  discussion  of  the  method 
of  images,  the  procedures  developed  in  this  section  are  limited  to  cases  in  which  the  incident 
radiation  propagates  normal  to  the  plane  of  the  mirror.  This  limitation  is  not  fundamental  to  the 
method  and  will  be  removed  in  a  future  work.  Section  III  reviews  two  simple  approximation 
techniques  that  have  been  used  in  previous  studies  of  light  scattering  from  a  particle  on  a  minor. 
These  methods  are  based  on  Mie  theory.  It  is  important  to  test  the  accuracy  of  these  approximations 
by  comparing  them  to  exact  results.  In  section  IV,  we  solve  several  test  cases  and  compare  the 
exact  results  to  the  results  obtained  by  the  two  approximation  procedures  discussed  in  section  JH. 
Section  V  ends  the  report  with  a  few  concluding  remarks. 
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II.  THEORY 


Consider  the  system  shown  in  Fig. la.  This  system  is  more  general  than  the  special  case  to 
be  considered  in  this  report.  We  begin  with  this  system  because  the  method  of  images  holds  for  the 
general  case  and  we  did  not  want  to  limit  the  discussion  to  a  specialized  case. 

The  z  a  0  plane  is  assumed  to  be  the  surface  of  a  perfect  conductor  (infinite  conductivity). 
One  or  more  particles  are  located  on  or  above  this  plane  in  the  region  z  >  0.  The  particles  can  have 
arbitrary  shapes  and  nonuniform  composition.  The  electrical  characteristics  of  the  system  are 
completely  described  by  specifying  the  dielectric  constant  e(r),  the  conductivity  a(r),  and  the 
magnetic  permeability  ji(r)  as  functions  of  the  position  coordinate  r  throughout  the  region  z  >  0.  In 
the  region  of  space  not  occupied  by  the  particles,  it  is  assumed  that  e(r)  =  1,  p(r)  =  1,  and  c(r)  =  0 
(Gaussian  units).  Inside  the  particle,  these  functions  take  on  the  values  characteristic  of  the 
particle.  On  the  z  =  0  plane,  the  conductivity  is  assumed  to  be  infinite.  A  plane  wave,  E^,  is 
incident  on  the  system.  This  gives  rise  to  a  reflected  plane  wave,  Eref,  and  a  scattered  wave,  E^. 
In  the  region  outside  of  the  particles,  the  total  electric  field  is  the  sum  of  these  three  contributions, 

E  =  Einc  +  Eref+Escat- 

The  electromagnetic  field  has  a  harmonic  time  dependence  given  by  exp(-ia>t).  In  the 
region  above  the  reflecting  plane,  it  must  satisfy  the  time  harmonic  form  of  Maxwell’s  equations 


VxH+ik(e  +  i4sa)E  =0 

V  x  E  -  ikpH  =  0 

V  •  (e  E)  =  0 

V  .(|iH)*0 


(1) 


where  k  =  co/c,  with  c  equal  to  the  velocity  of  light  The  field  must  also  satisfy  the  following 
boundary  conditions  on  the  z  =  0  plane: 


n  x  E  =  0 
n  H  =  0 


(2) 
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MIRROR  SURFACE 
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Fig.  1.  Illustration  of  the  method  of  images,  (a)  Physical  system, 
(b)  Equivalent  "method  of  images"  system. 


where  n  is  the  unit  vector  normal  to  the  plane.  These  conditions,  simply  stated,  say  that  there  can 
be  no  tangential  component  to  the  electric  field  and  no  normal  component  to  the  magnetic  field  at 
z  =  0. 

The  solution  to  this  problem,  stated  in  its  present  form,  is  quite  difficult  to  solve.  A 
technique  very  similar  to  the  method  of  images  in  electrostatic  problems  can  be  used  to  reformulate 
the  problem  in  a  much  simpler  form  without  explicit  boundary  conditions  at  the  surface  of  the 
mirror. 


A.  METHOD  OF  IMAGES 

The  method  of  images  is  intuitively  obvious.  An  observer  looking  down  on  the  system 
will  see  the  particles  and  mirror  images  of  the  particles  and  also  the  source  of  the  incident  wave 
(e.g.,  a  laser)  and  a  mirror  image  of  the  source.  The  image  of  the  source  is  the  apparent  origin  of 
the  reflected  wave.  This  apparent  picture  is  turned  into  reality  by  creating  the  system  shown  in 
Fig.  lb,  where  the  mirror  has  been  removed,  the  images  of  the  particles  have  been  replaced  by  real 
particles,  and  the  image  of  the  source  has  been  replaced  by  a  real  source.  The  image  source  is  the 
origin  of  the  plane  wave  E^.  By  definition,  Ej„c  =  Eref  in  the  region  z  >  0.  The  functions  e,  a, 
and  p,  in  this  model,  are  symmetric  with  respect  to  the  z  coordinate;  i.e.,  e(x,y,z)  =  e(x,y,-z), 
o(x,y,z)  =  o(x,y,-z)  and  p(x,y,z)  =  p(x,y,-z). 

By  the  principle  of  linear  independence,  each  of  the  incident  waves,  Einc  and  E,nc,  can  be 
considered  separately.  Einc  gives  rise  to  the  scattered  wave  Escat,  and  Einc  gives  rise  to  Escaf  in 
each  of  these  separate  scattering  problems,  a  plane  wave  is  scattered  by  a  muitiparticle  system.  It  is 
obviously  not  necessary  to  solve  both  of  these  problems  since  the  solutions  are  related  by 
symmetry.  Assume  that  we  have  solved  the  scattering  problem  for  the  field  E  =  Einc  +  Escat .  Write 
the  solution  in  Cartesian  coordinates 

E(x,y,z)  =  Ex(x,y,z)ex  +  Ey(x,y,z)ey  +  Ez(x,y,z)ez  (3a) 

where  ex,  ey,  and  ez  are  unit  vectors  along  the  x,  y,  and  z  axes,  respectively.  The  accompanying 
magnetic  field  is  written  similarly 

H(x,y,z)  =  Hx(x,y,z)ex  +  Hy(x,y,z)ey  +  Hz(x,y,z)ez  (3b) 
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These  fields  are  a  solution  to  the  scattering  problem  and  therefore  must  satisfy  Maxwell's 
equations.  We  now  claim  that  the  field  E  =  +  Escat  and  its  accompanying  magnetic  field  H  are 

related  to  E  and  H  by  the  expressions 

E(x,y,z)  =  -Ex(x,y,-z)ex  -  Ey(x,y,-z)ey  +  Ez(x,y,-z)e2 

(4) 

H(x,y,z)  =  Hx(x,y,-z)ex  +  Hy(x,y,~z)ey  -  Hz(x,y,-z)ez 


This  claim  is  justified  by  the  following  argument  Insert  the  above  expressions  for  E  and  H  into 
Maxwell's  equations  and  use  the  symmetry  of  the  functions  e(x,y,z),  a(x,y,z),  and  p(x,y,z)  with 
respect  to  z.  From  this  analysis,  we  see  that  if  E(x,y,z)  and  H(x,y,z)  satisfy  Maxwell’s  equations, 
then  E(x,y,z)  and  H(x,y,z)  will  also  satisfy  these  equations.  The  total  electric  and  magnetic 
fields  for  the  scattering  problem  shown  in  Fig.  1  b  are  obtained  by  summing  the  two  linearly 
independent  contributions 


Etoui  -  E  +  E 


H  I0tal  =  H  +  H 


(5) 


It  is  obvious  from  Eqs.  (3)  and  (4)  that  the  total  fields  Ewtal  and  Hwul  satisfy  the 
boundary  conditions  given  by  Eq.  (2).  They  also  satisfy  Maxwell's  equations.  Therefore,  in  the 
region  z  >  0  these  fields  are  a  solution  to  the  scattering  problem  shown  in  Fig.  la. 


B.  MULTIPOLE  EXPANSION  SOLUTION 

The  physical  problem  to  be  considered  in  this  report  is  illustrated  in  Fig.  2.  The  center  of  a 
spherical  particle  of  radius  R  is  located  a  distance  d  (where  d  >  R)  above  a  perfectly  conducting 
plane  surface.  The  center  is  assumed  to  lie  on  the  z  axis  and  the  conducting  surface  to  coincide  with 
the  z  =  0  plane.  A  circularly  polarized  plane  wave  traveling  in  the  negative  z  direction, 

Ejnc  =  -(ex  +  iey)  e  * kz,  is  incident  on  this  system.  By  the  method  of  images,  this  physical 
problem  is  converted  to  the  two-particle  scattering  problem  shown  in  Fig.  3.  In  addition  to  the 
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-(ex  +  iey)e-'kz 
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images,  to  the  problem  shown  in  Fig.  2. 


incident  wave,  an  image  wave,  traveling  in  the  positive  z  direction,  Emc  =  (ex  +  i  ey)  e1  kz,  is  also 
present.  The  total  incident  wave  for  the  two-particle  problem  is  the  sum  of  these  two  waves 

Ei=  (ex  +  iey)  (e^-e-i  (6) 

This  wave  satisfies  the  boundary  conditions  given  by  Eq.  (2)  on  the  z  =  0  plane. 

In  the  region  outside  the  spheres,  the  total  electric  field  is  conveniently  written  as  the  sum 
of  the  incident  wave  and  two  contributions  to  the  scattered  wave 

E=Ei  +  E(s1,  +  E®  (7) 

The  two  scattered  waves,  E* }  and  E(s2) ,  propagate  radially  outward  from  the  centers  of  die  two 
spheres,  Oj  and  02 ,  respectively  (see  Fig.  3). 

The  electric  fields  can  be  expanded  in  terms  of  the  vector  spherical  wave  functions 


M$m  =  zg>(kr)eM>xn,m(Q) 

and  (8) 

N&  -  4 '(ta-)]  Y„.m<8)  +  zS>(kr)  Zm(8)j 

where  zjJ)  is  a  spherical  Bessel  function  of  type  jn,h^  for  j  =  1,3  respectively.  The  vector 
functions  X,  Y,  and  Z  are  defined  by 


Xn,m(9)  —  i  ^n,m(®)  e0  '  e$ 


Yn,m(0)  —  ^n,m(0)  e0  +  i  ftn,m(0) 


(9) 


Zn,m(0)  =  n(n+l)  PjJ1  (cos8)er 


where 

^n.m(9)  =  ~m—  P?(cos  8) 
sin(O) 

(10) 

Vm(0)=— P?(cos  8) 

ae 
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The  function  P?(x)  is  the  associated  Legendre  polynomial,  and  (er,  ee,  e$)  are  unit  orthogonal 
vectors  associated  with  the  spherical  coordinates  r,8,<j). 

In  the  analysis  that  follows,  only  the  functions  Mn  m  and  Nnjn  with  m  =  1  will  be  needed. 
Therefore,  to  simplify  the  notation,  the  subscript  m  =  1  will  not  be  explicitly  written.  It  wall  be 
understood  that  Mn  =  Mn  l  and  Nn  =  Nnl.  Using  the  standard  formulas  for  the  expansion  of  a 
linearly  polarized  plane  wave,10'11  we  can  derive  the  following  formula  for  the  expansion  of  a 
circularly  polarized  plane  wave 


oo 

(ex  +  i  ey)  e±  * k z  =  +£  (±i)n+1  ^^[M<nl>(r,0,<|))  ±  N^(r,0,0)]  (11) 


In  the  discussions  that  follow,  the  coordinates  r,0,4>  are  defined  with  respect  to  coordinate 
axes  with  origin  at  0  (see  Fig.  3).  A  point  in  this  coordinate  space  will  also  be  represented  by  the 
boldfaced  position  vector  r.  In  addition,  it  will  be  useful  to  introduce  two  other  coordinate  systems 
with  origins,  0X  and  02,  located  it  the  centers  of  the  two  spheres.  The  coordinates  in  these  systems 
will  be  identified  with  the  subscripts  1  and  2,  respectively  [i.e.,  rj  or  (rp  0j,4>j );  j  =  1,2].  These 
three  sets  of  coordinates  are  related  to  each  other  by  translations  along  the  z  axis. 

The  incident  wave,  defined  by  Eq.  (6),  can  be  expanded  in  terms  of  spherical  wave 
functions  centered  at  0! 


Ei(rj)=  X  [PnM^rO  +  qnN^rO]  (12) 

n  =  l 


The  expansion  coefficients  in  this  equation  can  be  obtained  with  the  aid  of  Eq.  (11).  They  are 


pn  =  -in+1  2n  ±--1-  [e-*d  -  (-1)"  e*d] 
n(n+l) 


qn  =  -in+l  [e-ikd  +  (-[)"  ei*d] 


(13) 


where  d  is  the  displacement  of  the  center  of  the  sphere  from  origin  0  (see  Fig.  3).  In  addition,  the 
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scattered  waves  EjP  and 


can  be  expanded  around  centers  Oj  and  O2,  respectively. 


Ep)(rj)  =  X  [ag>  M^(rj)  +  bp>  Njftrj)] ;  j  =  1,2  (14) 

n  = 


The  scattering  coefficients  a^\  a£\  and  t&2)  are  to  be  determined  by  the  calculation.  This  task 

is  made  somewhat  easier  by  exploiting  the  mirror  symmetry  of  the  system  and  the  boundary 
conditions  on  the  z  =  0  plane.  These  conditions  imply  the  following  relations  between  the 
scattering  coefficients  (see  appendix  A): 


41)  =  -(-Dn42> 

b£>=  C-iy'bP 


(15) 


To  carry  out  the  calculations  of  the  scattering  coefficients,  it  is  convenient  to  expand  the 
fields  in  the  coordinate  system  with  origin  at  Oj.  The  fields  E;  and  El1)  are  already  in  this  form. 
El2)  is  not,  but  can  be  converted  to  this  form  by  means  of  the  translation-addition  theorem  for 
vector  spherical  wave  functions.  A  general  formulation  and  discussion  of  this  theorem  is  given  in 
references  8  and  9.  The  special  form  of  this  theorem,  needed  for  the  present  problem,  is  the 
following: 


Ml3>(r2)=  £  [An.n-  + 

ri  =  1 

(16) 

00 

Nj?>(r2)=  £  [Ahj,*  N^.J(n)  +  Bn,n*  M^Vl)] 

ri  SK  1 

The  formulas  for  the  expansion  coefficients  A^-  and  B,^  and  a  discussion  of  the  method  for 
calculating  these  quantities  are  presented  in  appendix  B. 

The  electric  field  in  the  region  outside  the  spheres  is  given  by  Eq.  (7).  The  multipole 
expansion  of  this  field  around  origin  Oj  can  be  derived  with  the  aid  of  Eqs.  (12),  (14),  and  (16). 
The  result  is 
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E  =  X  {  [Pn  +  X  (A„.na<?  +  B„-,nb(n2))]Mj11)(r1)  +  a^M^n) 

n  n' 

(17) 

+  [qn  +  Z  (B„-,„a<J>  +  An.)nb(n2)>]Nj11)(ri)  +  btf>  Nj?>(n) } 

n* 


This  equation  is  of  the  form 

E  =  X(an  M£>  +  a^Mk3>  +  pn  N11*  +  b^^N^)  (18) 

n 


The  ratios  of  the  amplitudes  of  the  scattered  wave  multipoles  (M(3)  and  N(3))  to  the  incident  wave 
multipoles  (M(1)  and  N{1))  are 


u 


Ctn 


(19) 


The  quantities  and  vn  are  elements  of  the  so-called  scattering  T  matrix12  for  the  particle. 
For  a  homogeneous  sphere,  and  vn  are  the  Mie  theory  coefficients  for  the  TE  and  TM  scattering 
modes,  respectively.  The  formulas  for  u„  and  vn,  for  the  homogeneous  dielectric  sphere,  are  given 
in  appendix  C.  [  It  should  be  noted  that  it  is  not  necessary  for  the  particle  to  be  homogeneous.  If 
the  index  of  refraction  is  a  function  of  the  radius,  u„  and  vn  can  still  be  calculated  by  an  appropriate 
theory.  For  example,  the  theory  of  Aden  and  Kerker 13  or  Bhandari14  can  be  used  for  a  layered 
sphere,  or  the  theory  of  Wyatt15  can  be  used  for  a  sphere  with  continuous  radially  varying  index  of 
refraction.] 

Assume  now  that  the  Mie  theory  coefficients,  Uj,  and  vn ,  have  been  calculated  for  the 
particle.  Substitute  the  expressions  for  otj,  and  pn ,  obtained  from  Eq.  (17),  into  Eq.  (19)  to  obtain 


a£l)  =  un 


bff>»vn 


(20) 
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The  sums  over  n’  in  Eq.  (20)  converge  and  can  be  truncated  to  the  range  1  <  ri  <  N,  where  N  is 
the  number  of  modes  required  for  convergence  of  a  Mie  theory  calculation  for  scattering  from  the 
isolated  spherical  particle. 

These  equations  can  be  cast  in  a  much  more  elegant  and  easy-to-use  form  by  expressing  all 
the  quantities  as  either  square  matrices  or  column  vectors.  The  quantities  a^\  b^,  a^2),  b^2),  pn, 
and  are  elements  of  the  (length  N)  column  vectors:  aO),  b(1),  a(2\  b<2),  p,  and  q.  The 
quantities  A,,^  and  are  elements  of  the  N  by  N  square  matrices  A  and  B.  The  scattering 
coefficients  %  and  vn  are  elements  of  the  diagonal  matrices  u  =  [Uj,  and  v  =  [vn  It  is 
also  useful  to  define  the  diagonal  matrix  g  =  [(-l)n  5^].  The  set  of  2N  linear  equations,  Eq.  (20), 
can  now  be  written  in  the  super  matrix  form 


/  AT  +  gu1  BT  \/a<2>\  /P) 

\  Bt  At  -  gv1  /l  b(2)  /  Q 


(21) 


where  we  have  taken  account  of  Eq.  (15)  to  eliminate  a0>  =  -  ga(2>  and  bO)  =  gb(2>  from  the 
equations.  The  superscript  T  indicates  the  transpose  matrix.  Standard  linear  equation-solving 
subroutines  are  available  on  most  computers  to  solve  this  matrix  equation  for  the  scattering 
coefficients  a<2>  and  b(2>. 


C.  CROSS  SECTIONS 

The  scattered  wave  is  the  sum  of  the  two  scattering  components 


Escat  =  Ei1)  +  Ei2>  (22) 

The  asymptotic  form  of  this  wave  can  be  evaluated  in  the  limit  r  — >  <»  with  the  aid  of  Eq.  (14)  and 
the  following  asymptotic  formulas  for  the  spherical  wave  functions: 


MV=(-ir^xnm« 

N^  =  i(-i)n^Yn(e)ei«!> 


(23) 
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The  result  is 


Escat  =  f^  £  (-iytaOJXnCeo  +  ibtfJYnCeo] 

1  n  =1 

(24) 

+  ^  ei02  £  (-i)n  [42)  Xn(02)  +  i  *42)  Yn(02>] 

^2  n=l 


This  formula  is  expressed  in  terms  of  the  two  coordinates  systems,  (rlt  0j,  (J^)  and 
(r2, 02,  <(>2)-  To  be  useful,  it  must  be  expressed  in  terms  of  a  single  set  of  coordinates  (r,0,  <>).  The 
asymptotic  transformation  relations  between  the  coordinates,  in  the  limit  r  — >  «>,  are 


ri  =  r  +  d  cos(0) 
r2  =  r  *  d  cos(0) 

01  =02  =  0  (25) 

4>1  =  <J>2  *  <}> 

Substitute  these  transformation  relations  into  Eq.  (24)  to  obtain  the  scattered  wave  in  terms  of  r,  0 
and  <(>.  This  wave  can  be  written  in  the  form 

Escar  =  Sfc  e»  [Se(9)  e9  +  i  S0(0)  e J  (26) 


where  the  vector  components  of  the  scattering  amplitude  are  given  by 

Se(0)  =  c^dcos(Q)  S^(0)  +  e-^«>s(0)  Sq2)(0) 

S^(0)  =  eikdcos(0>  S^!)(0)  +  e-ikdcos<0)  S®(8) 

and  where 

N 

S ®(0)  =  -  £  (-i)n+1  [aS>  Kn  (0)  +  bSP  t„  (0)] 

n  =  1 

S®(0)  =  -  £  (-i)n+1  [a?  tn  (0)  +  bfi>  jc«  (0)j 

n  =  1 

for  j  =  1,2. 


(27) 


(28) 
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The  formulas  for  the  differential  scattering  cross  section  can  be  taken  directly  from  the 
similar  formulas  obtained  in  Mie  theory.10-11  If  the  incident  beam  is  right  or  left  circularly 
polarized  or  if  it  is  unpolarized,  the  differential  scattering  cross  section  is  given  by 


a(0)  =  ^[lSe(0)|2  +  lS*(8)|2] 


(29) 


If  the  incident  beam  is  plane  polarized,  with  the  electric  field  vector  parallel  to  the  x  axis,  the 
differential  cross  section  is  given  by 

0(0)  =  \  [I  S0  (0)  i  2cos2(q)  + 1 S*  (0)  1 2sin2(<j>)]  (30) 


If  the  electric  field  is  polarized  parallel  to  the  y  axis,  the  factors  cos2(<$>)  and  sin2^)  in  the  above 
formula  are  interchanged.  The  angle  0  in  these  formulas  is  restricted  to  the  region  above  the 
reflecting  plane,  i.e.,  0  <  8  <  it/2. 

The  total  scattering  cross  section  is  the  integral  of  the  differential  cross  section  over  the  2tc 
steradian  solid  angle  above  the  z  =  0  plane.  For  either  of  the  differential  cross  section  formulas, 
given  in  Eqs.  (29)  or  (30),  the  result  is 

JfJt/2 

[l  Se(0)  P  + 1  S^(0)  |2]  sin(0)  d8  (31) 

0 

In  Mie  theory,  it  is  possible  to  substitute  the  series  expansions  for  Se(0)  and  S^(0)  into  Eq. 
(31)  and  then  use  the  orthogonality  properties  of  the  functions  7^(0)  and  xn(8)  to  obtain  a  simple 
formula  for  the  total  scattering  cross  section.  This  approach  is  not  feasible  in  the  present  case 
because  the  functions  Se(0)  and  S^O)  are  much  more  complicated  than  in  Mie  theory.  Instead,  the 

scattering  cross  section  is  calculated  using  the  general  relation 

Cscat  =  Cext  *  Cabs  (32) 

where  is  the  extinction  cross  section  and  is  the  absorption  cross  section. 

The  extinction  cross  section  can  be  computed  with  the  aid  of  the  optical  theorem11 


where  S(0)  =  Se(0)  =  S$(0)  is  the  scattering  amplitude  in  the  forward  direction  6  =  0.  Use  the 
relations  7t„(0)  =  tn(0)  =  n(n+l)/2  and  Eqs.  (27)  and  (28)  to  obtain 

=  X  (-i)"n(n+l)[e‘kd(a(n1)  +  b(Ql))  +  e-^a^  +  b(n2))  ]  (34) 

k2  n=  l 

The  absorption  cross  section  is  calculated  by  utilizing  formulas  for  power  absorption 
developed  in  Mie  theory.  The  Mie  theory  formula  for  the  power  absorbed  by  spherical  particle 
scattering  the  plane  wave  (ex  +  i  ey)  e**2  is 

P  =  F X  (2n  +  1)  [(]un|2  +  Re  un)  +  (|  vn|2  +  Re  vn)]|  (35) 

\  kz  n -  1  | 

where  F  is  the  energy  flux  of  the  plane  wave  and  the  quantity  in  braces  is  the  absorption  cross 
section.  We  now  interpret  this  formula  in  terms  of  the  contributions  of  each  of  the  multipole 
components  of  the  plane  wave  to  the  total  power  absorbed  by  the  particle.  Define  P(jM>  and  P^  to 
be  the  power  absorbed  from  the  n'th  TE  and  TM  multipole  components,  respectively,  of  the 
incident  plane  wave 


P^>  =  ^ (2n+  x>  (lunl2  +  Re  uj 
P^  =  ij-  ^  (2n+  i)  (1  vn|2  +  Re  vn) 


(36) 


The  expansion  of  the  plane  wave  is  given  by  Eq.  (1 1).  The  amplitude  of  each  of  the 
multipole  components  of  this  wave  is 


fn  =  - 


jn+i  2n-f-l 

n(n+l) 


(37) 


The  power  absorbed  from  each  multipole  component  is  proportional  to  the  square  of  the  amplitude 
fn .  Therefore,  to  obtain  the  power  absorbed  from  a  multipole  of  unit  amplitude,  we  must  divide 
the  results  given  in  Eq.  (36)  by  lfn  I2.  The  result  is 
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(38) 


n(n+l) 


2 


L2n+  U 


PkM) 


n(n+l) 
2n  +  1 . 


2 


The  wave,  which  is  incident  on  particle  1,  has  a  multipole  expansion  that  can  be  deduced  from  Eq. 
(18): 

E-£fa.Mg>  +pnNi»)  (39) 


The  expansion  coefficients,  obtained  from  Eq.  (19),  are 


a 


u„ 


(40) 


The  total  energy  absorbed  by  particle  1  is  obtained  by  multiplying  P„M)  and  [defined  in  Eq. 

(38)]  by  the  squares  of  the  amplitudes  |  an  |2  and  j  (3n  f  and  summing  over  all  the  terms.  The 
absorption  cross  section  is  then  obtained  by  dividing  this  result  by  the  energy  flux  F.  The  final 
result  is 


N 


c  K  _  _  2tl  V  [n(n+l)F 
°  k2  i  2n  + 1 

*■  n  =  1 


4^-  (hni2  +  Re  un)  +  l^l  (| vnp  +  Re  vn) 


(41) 


These  formulas  were  tested  by  comparing  the  scattering  cross  section  obtained  by 
numerically  integrating  the  differential  cross  section,  given  by  Eq.  (31),  with  the  results  obtained 
using  Eqs.  (32),  (34),  and  (41).  The  tests  were  carried  out  using  both  real  and  complex  values  for 
the  index  of  refraction.  The  results  were  in  complete  agreement,  within  the  numerical  accuracy  of 
the  calculation. 
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III.  APPROXIMATION  METHODS 


Approximation  methods  have  been  used  in  previous  studies  to  calculate  the  differential 
scattering  cross  section  for  a  particle  on  a  conducting  plane.2*6  Two  of  these  approximations  will 
be  reviewed  in  this  section.  In  addition,  a  technique  known  as  the  order  of  scattering  (OS)  method 
will  also  be  presented,  and  the  relationship  of  this  method  to  one  of  the  approximations  will  be 
discussed. 


A,  FORWARD  SCATTER  MIE  APPROXIMATION 

This  method  is  illustrated  in  Fig.  4.  Mie  theoiy  is  used  to  calculate  the  scattering  from  the 
spherical  particle.  The  forward  scattered  rays  are  reflected  from  the  mirror  back  to  the  detector.  The 
back  scattered  rays  are  usually  much  weaker  than  the  forward  scattered  rays  and  are  therefore 
neglected.  The  mirror  is  assumed  to  have  no  effect  on  the  Mie  scattering  cross  section.  The  only 
effect  of  the  mirror  is  to  reflect  the  forward  scattered  rays.  Young2  developed  this  approximation 
and  obtained  good  agreement  with  experimental  results.  Nahm  and  Wolfe6  also  used  this  method, 
which  they  refer  to  as  the  unobstructed  reflection  model. 


B  .  SINGLE  SCATTER  APPROXIMATION 

The  single  scatter  approximation  (SSA)  method  is  also  based  on  Mie  theoiy.  This  method 
takes  account  of  all  of  the  ways  in  which  a  ray  in  the  incident  beam  can  interact  only  once  with  the 
particle  and  be  deflected  to  the  detector.  The  four  ways  this  can  happen  are  illustrated  in  Fig.  5. 
This  method  neglects  all  processes  in  which  a  ray  is  scattered  two  or  three  or  more  times  by  the 
particle.  This  would  happen,  for  example,  if  a  ray  were  deflected  by  the  particle  toward  the  mirror, 
then  were  reflected  back  and  struck  the  particle  a  second  time,  and  then  were  deflected  to  the 
detector.  Such  multiple  scattering  processes  are  ignored  in  this  approximation.  The  amplitudes  of 
the  four  scattered  rays  (waves)  shown  in  Fig.  5  are  added,  taking  account  of  the  phase  differences 
due  to  the  different  distances  that  the  incident  rays  must  travel  to  reach  the  particle,  the  different 
distances  that  the  scattered  rays  must  travel  from  the  particle  to  the  detector,  and  the  phase  shift  of 
it  radians  caused  by  reflection  from  the  mirror.  Only  a  single  Mie  theory  calculation  is  required  to 
carry  out  this  calculation;  therefore,  the  method  is  very  efficient.  This  method  is  almost  identical  to 
the  double  interaction  model  (DIM)  of  Nahm  and  Wolfe.6  The  DIM  includes  an  "area  factor,” 
which  is  not  included  in  the  SSA  model. 
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Fig.  4.  Forward  scatter  Mie  approximation. 
Backscatter  is  ignored. 


Fig.  5.  Single  scatter  approximation  (SSA).  The  amplitudes  of  these 
four  single  scatter  events  are  added,  taking  account  of  phase 
differences  due  to  the  different  path  lengths. 
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The  SS  A  can  be  considered  to  be  the  first  term  of  an  infinite  series  expansion  in  which  the 
second  term  takes  account  of  all  double  scattering  events,  the  third  term  takes  account  of  all  triple 
scattering  events  and  so  on.  The  general  technique  of  summing  these  multiple  scattering  events  is 
known  as  the  OS  method.  This  method  and  its  relationship  to  the  exact  theory  will  be  discussed 
next. 


C.  ORDER  OF  SCATTERING  METHOD 

Most  of  the  computational  effort  in  solving  the  exact  scattering  problem  is  due  to  Eq.  (21). 
This  matrix  equation  is  of  the  form 


(u  1  -  c)f  =  q 


where 


(42) 


and 


The  formal  solution  to  Eq.  (42)  is 

F  =  (U‘1  -  C)-1  Q 


(43) 


(44) 


(45) 


(46) 


(47) 


The  inverse  matrix  in  this  equation  can  be  evaluated  by  the  following  infinite  series  expansion 
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(u*1-  c)"1  =  u  +  ucu  +  ucucu  + 


(48) 


If  this  series  converges  rapidly,  it  is  a  useful  procedure  for  solving  Eq.  (42).  This  technique  was 
developed  by  Fuller  and  Kattawer16  and  is  used  to  calculate  scattering  from  an  ensemble  of 
spheres.  The  first  term  of  the  series  represents  the  single  scattering  contribution,  the  second  term  is 
the  double  scattering  contribution,  and  so  on.  If  we  retain  only  the  first  term,  the  result  is  the  SSA 
method.  The  scattering  coefficients  for  this  case  can  be  easily  calculated: 


41)  =  unPn 


b(i1)  =  vn  qn 


(49) 


This  requires  a  single  Mie  theory  calculation  to  obtain  u„  and  vn. 

Based  on  the  simple  ray  optics  picture  of  the  scattering  process,  one  would  expect  the 
probability  of  double  scattering  events  to  be  small  relative  to  single  scattering  events  whenever  the 
scattering  cross  section  is  small  (i.e.,  for  small  particles)  or  when  the  distance  d  of  the  panicle  from 
the  mirror  is  large.  Under  these  circumstances,  the  single  scatter  approximation  is  expected  to  be 
accurate  and  should  converge  to  the  exact  result  in  the  limit  that  the  particle  size  approaches  zero  or 
the  distance  d  approaches  infinity.  In  the  next  section,  we  will  see  that  these  expectations  are 
apparently  borne  out 
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IV.  RESULTS 


The  results  of  several  test  calculations,  which  compare  the  predictions  of  the  exact  theory 
with  the  predictions  of  the  approximation  methods  (the  forward  scatter  Mie  approximation  and  the 
single  scatter  approximation  method),  are  presented.  These  results  give  a  preliminary  estimate  of 
the  accuracy  of  the  approximation  methods.  The  cross  sections  are  given  in  reduced  units  of  area. 
The  reduced  unit  of  length,  in  this  system,  is  defined  to  be  the  wavelength  of  the  scattered  light 
(i.e.,  in  these  units,  X  =  1).  Results  for  particles  of  three  different  sizes  are  shown  in  Figs.  6, 7, 
and  8.  The  radii  of  these  particles,  in  reduced  units,  are  R  =  0.2, 1.0,  and  2.5.  The  particles  are 
resting  on  the  surface  of  the  mirror,  therefore,  the  distance  parameters  are  given  by  d  =  R.  The 
index  of  refraction  is  n  =  1.46.  As  anticipated,  the  SSA  results  are  fairly  good  for  the  small 
particle  size. 

In  each  of  the  cases  shown  in  Figs.  6-8,  the  Mie  theory  approximation  gives  cross  sections 
for  small  angle  scattering  (0  <  0  <  20°)  that  are  smaller  than  the  exact  cross  sections.  To  dispel 
any  notion  that  this  might  always  be  true,  we  show  Fig.  9,  where  the  the  index  of  refraction  has 
been  changed  to  n  =  1.3.  The  other  parameters  are  the  same  as  in  Fig.  7,  i.e.,  R  =  1.0  and  d  =  1.0. 
The  Mie  theory  results  are  now  much  larger  than  the  exact  results  for  low  angle  scattering. 

Figure  10  shows  results  for  a  particle  that  is  suspended  above  the  plane  of  the  mirror.  The 
parameters  for  this  case  are  R  =  1.0,  n  =  1.3,  and  d  =  5.0  (i.e.,  the  center  of  the  sphere  is  5.0 
units  above  the  mirror).  The  SSA  method,  as  expected,  is  a  fairly  good  approximation  for  this 
case. 
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CROSS  SECTION 


Fig.  6.  Cross  section  for  light  scattering  from  a  particle  on  a  mirror.  R  =  0.2,  d  =  0.2, 

n  =  1.46.  Comparison  of  exact  solution  with  the  forward  scatter  Mie  approximation 
(Mie)  and  the  single  scatter  approximation  (SSA). 
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CROSS  SECTION 


Fig.  7.  Cross  section  for  light  scattering  from  a  particle  on  a  mirror.  R  =  1.0,  d  =  1.0, 

n  =  1.46.  Comparison  of  exact  solution  with  the  forward  scatter  Mie  approximation 
(Mie)  and  the  single  scatter  apnroximation  (SSA). 
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CROSS  SECTION 


SCATTERING  ANGLE 


Fig.  8.  Cross  section  for  light  scattering  from  a  particle  on  a  mirror.  R  =  2.5,  d  =  2.5, 

n  =  1 .46.  Comparison  of  exact  solution  with  the  forward  scatter  Mie  approximation 
(Mie)  and  the  single  scatter  approximation  (SSA). 
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Fig.  9.  Cross  section  for  light  scattering  from  a  particle  on  a  mirror.  R  =  1.0,  d  =  1.0, 

n  ~  1-3*  Comparison  of  exact  solution  with  the  forward  scatter  Mie  approximation 
(Mie)  and  the  single  scatter  approximation  (SSA). 
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Fig.  10.  Cross  section  for  light  scattering  from  a  particle  near  a  mirror.  R  =  1.0,  d  =  5.0, 

n  =  1.3.  Comparison  of  exact  solution  with  the  forward  scatter  Mie  approximation 
(Mie)  and  the  single  scatter  approximation  (SSA). 
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V.  CONCLUDING  REMARKS 


This  report  describes  a  practical,  efficient  computational  method  for  calculating  light 
scattering  from  a  spherical  particle  on  a  conducting  surface.  The  computational  efficiency  of  this 
method  makes  it  possible  to  eliminate,  for  the  most  part,  the  need  to  use  approximation  methods. 
We  are  presently  using  this  procedure  to  calculate  mirror  Bidirectional  Reflectance  Distribution 
Functions  (BRDF)2  due  to  various  distributions  of  particulate  contamination  on  the  surface  of  a 
mirror. 

The  present  form  of  the  theory  is  restricted  to  the  case  in  which  the  incident  light  propagates 
in  a  direction  normal  to  the  surface  of  the  mirror.  This  limitation  will  be  removed  in  future  work. 

The  listing  of  a  Fortran  computer  code  that  implements  the  procedures  outlined  in  this 
report  is  presented  in  Appendix  D  . 
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APPENDIX  A 

SCATTERING  COEFFICIENT  SYMMETRY  RELATIONS 


The  symmetry  relations  expressed  by  Eq.  (15)  in  the  text  are  a  consequence  of  the 
symmetry  of  the  system  and  the  boundary  conditions  on  the  z  =  0  plane.  The  boundary  conditions 
imply  that  at  any  point  on  this  plane,  the  component  of  the  electric  field  that  is  parallel  to  the  plane 
must  vanish.  This  applies  to  the  total  field  given  by  Eq.  (7).  Since  the  incident  field,  Ei;  satisfies 
this  condition  separately,  it  follows  that  the  scattered  field,  E^  =  E(s!)  +  E$?,  must  also  satisfy  the 
condition  separately. 

To  carry  out  the  analysis,  we  choose  an  arbitrary  point  P  on  the  z  =  0  plane.  Assume  that 
this  point  has  coordinates  r^Qj,^  with  respect  to  the  coordinate  system  centered  at  Oj.  Then,  the 
coordinates  of  P  with  respect  to  the  02  coordinate  system  are  given  by 


r2  =  rj 

02  =  n  -  6]  (Al) 

02  =  01 

At  point  P,  the  unit  vectors  e$,  e<^  and  e^2  are  collinear  and  lie  in  the  z  =  0  plane.  Thus, 
the  e^,  component  of  the  total  scattered  field  is  obtained  by  algebraically  adding  the  and  e$2 
components  of  the  E^^  and  E^2)  fields.  The  e$,  and  e^2  components  of  the  scattered  field  can  be 
obtained  with  aid  of  Eqs.  (8),  (9),  and  (14).  The  result  is 

l$i\  =  -  X  exP(i0j)  *n(0j)]  -  i  [rj  zn  )(kfj)]|  [  HP  7tn(0j)]J  (A2) 

where  j  =  1,2.  The  boundary  conditions  at  P  require  that 

[EH, +[#>]*,-<>  (A3) 

This  relation  can  only  be  true  if  it  holds  for  each  individual  multipole  component  of  the  field.  This 
implies  that 


a^Tntfi)  +  42>xn(02)  =  0 

and  (A4) 

b^1>7tn(01)  +  b^)jTn(02)  =  O 
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Use  Eq.  (Al)  and  the  following  symmetry  relations  for  the  and  xn  functions 


Kn(TC-e)  =  -(-l)n7Cn(0) 
t„0c-0)  =  (-  l)nx„(8) 
to  obtain  the  desired  results  given  by  Eq.  (15)  in  the  text. 


(A5) 


34 


APPENDIX  B 

THE  TRANSLATION-ADDITION  THEOREM 

The  special  form  of  the  translation-addition  theorem  used  is  this  report  is  given  by  Eq.  (16) 
in  the  text.  This  formula  expresses  the  vector  spherical  wavefunctions  M(n3>(r2)  and  N(n3)(r2)  about 
origin  02  in  terms  of  wave  functions  and  N^Cri)  about  origin  0^  where  02  is  displaced  a 

distance  6  =  2d  along  the  z  axis  from  0j  (see  Fig.  3).  The  expansion  coefficients,  A„  n.  and  Bn  n., 
which  appear  in  Eq.  (16),  are  obtained  from  the  general  formulas  given  in  the  appendix  of 
reference  6: 

A„  n'  =  -in'-V¥#7rZ  i_V  Mn+1>  +  n’<n'+1)  -  v(v+1>]  a(n>n"v>  hv1}(k5)  (Bl) 
2n(n+l)rr 


B.  _  :n'-n 
n,n  -  1 


2n'+l 

2n'(n'+l) 


i-v 


(2ik5)a(n,n';v)ht1>(k5) 


(B2) 


where  h^(k8)  is  the  spherical  Hankel  function  of  the  first  kind  and  a(n,n';  v)  is  special  form  of  the 
Guant  coefficient17,  defined  by  the  following  product  of  two  3-j  symbols 


a(n,n';v)  =  (2v+l) 


n(n+l)  W2  /„  n-  v 


Ln’(n’+1) 


(n  n'  vl(n  n’  v\ 

\o  o  oMi  -1  0 1 


(B3) 


The  definition  of  the  3-j  symbol  involves  the  summation  of  many  factorial  terms.  As  a 
result,  a  straightforward  evaluation  of  a(n,n’;  v)  is  very  inefficient.  Bruning  and  Lo9  found  a  way 
to  overcome  this  difficulty  by  developing  a  recurrence  relation  for  a(n,n’;  v)  that  cycles  on  the 
index  v. 

We  did  not  use  this  method.  Instead,  we  found  it  more  convenient  to  use  a  procedure 
developed  by  Schulten  and  Gordon18  for  calculating  the  3-j  symbols.  This  procedure  is  also  based 
on  a  recurrence  relation  that  cycles  on  the  index  v.  It  was  more  convenient  because  we  could  make 
use  of  an  existing  computer  subroutine  available  to  us  to  do  the  calculations.  We  found  this 
approach  to  be  very  efficient. 
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APPENDIX  C 

FORMULAS  FOR  THE  uD  AND  vn  COEFFICIENTS 


The  Ujj  and  vn  coefficients  are  defined  by  Eq.  (19)  in  the  text  For  the  case  of  a 
homogeneous  sphere,  they  arc  identical  to  the  Mie  theory  "aj,  and  bn"  coefficients  for  the  TE  and 
TM  scattering  modes,  respectively.  Since  there  are  many  different  definitions  of  these  quantities  in 
the  literature,  it  seems  advisable  to  give  the  formulas  for  these  quantities  that  arc  consistent  with 
their  usage  in  this  report  The  definitions  we  use  are  the  same  as  those  defined  by  Stratton10 
(Section  9.25).  If  we  take  the  permeability,  ji,  of  the  particle  and  the  surrounding  medium  to  be 
the  same,  these  coefficients  are  given  by 


\j/n(mx)  y'n(x)  ~  m  yn(x)  ^(mx) 
Vn(mx)  5'n(x)  -  m£n(x)  v'n(mx) 

¥n(x)  V'n(mx)  -  myn(mx)  yn(x) 
£n(x)  V'n(mx)  -  m\j/„(mx)  £’n(x) 


(Cl) 

(C2) 


where  m  is  the  complex  index  of  refraction;  x  =  kR  is  the  size  parameter  of  the  particle,  the  prime 
means  to  take  the  derivative  with  respect  to  the  argument  of  the  function;  and  \yn(x)  and  ^(x)  are 
Riccati-Bessel  functions  defined  by  vyn(x)  =  x  jn(x)  and  ^n(x)  =x  h^fx). 
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APPENDIX  D 

FORTRAN  COMPUTER  CODE 


On  the  following  pages  is  a  listing  of  the  Fortran  computer  code  that  implements  the 
procedures  outlined  in  this  report.  The  program  runs  interactively  and  prompts  the  user  for  input 
data.  The  code  is  efficient  and  can  quite  reasonably  be  run  on  a  personal  computer.  The  output 
consists  of  the  differential  scattering  cross  section  as  well  as  the  extinction,  absorption,  and  total 
scattering  cross  sections  for  both  the  "spherical  particle  on  a  conducting  plane"  and  for  the 
"isolated  sphere."  The  scattering  results  for  the  "isolated  sphere"  (Mie  theory)  are  a  natural  by¬ 
product  of  the  calculation.  The  results  for  the  isolated  sphere  are,  by  definition,  the  "forward 
scatter  Mie  approximation"  cross  sections. 

Following  the  program  listing  is  a  sample  output  that  was  calculated  using  the  default  input 
parameters.  The  cross  sections  in  this  example  output  correspond  to  the  "exact"  (and  also  the 
"forward  scatter  Mie")  cross  sections  shown  in  Fig.  7.  (However,  many  more  points  were 
calculated  for  the  graph  than  are  shown  in  the  sample  output) 

The  code  is  self-contained  except  for  the  IMSL  Library  subroutine  LECT1C,  which 
appears  in  the  subroutine  FSOLVE.  LECTIO  is  a  complex  linear  equation  solver  that  solves 
Eq.  (21)  in  the  text  The  user  must  either  supply  this  subroutine  or  an  equivalent  subroutine. 
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PROGRAM  MIRMX ( i nput , output , tape  5=i nput , tape 6=output ) 
COMMON/RMIE/RAD, RFR, RFI, WL, KAY, LMAX 
COMMON/AB/  A  (60, 60) ,B (60, 60) 

COMMON/CMIE/U (100) ,V(100) 

CQMMON/GFMAT/G (120, 120) , F (120) 

COMMON/DIFF/NANG,THETA(400) ,DIFF (400) , DIFS (400) , QEXT, QSCA, QABS 
C0t®40N/CRS2/Q2EXT,  Q2ABS 
DIMENSION  TWRT(10) 

COMPLEX  A, B,U,V, G,F, CSOM, UA,  UB, VA, VB 
REAL  KAY 
C 

C  THIS  PROGRAM  CALCULATES  SCATTERING  FROM  A  SPHERICAL 
C  PARTICLE  ON  A  MIRROR.  IT  USES  THE  METHOD  OF  IMAGES  TO 
C  CONVERT  TO  AN  EQUIVALENT  2  PARTICLE  PROBLEM  AND  THEN  USES 
C  THE  METHOD  OF  BRUNING  AND  LO  TO  SOLVE  THIS  EQUIVALENT 
C  SCATTERING  PROBLEM.  (B.R.  JOHNSON,  MAY.  1990) 

C  MAJOR  CORRECTIONS  AND  REVISIONS  OCT.  1990) 

C  INTERACTIVE  INPUT  VERSION  DEC.  1990 
C 

C  PARAMETER  DEFINITIONS 
C  WL  =  WAVELENGTH 
C  RAD  =  PARTICLE  RADIUS 

C  DELANG  =  INTERVAL  BETWEEN  ANGLES  OF  DIFFERENTIAL  CROSS  SECTIONS. 

C  NANG  =  NUMBER  OF  DIFFERENTIAL  CROSS  SECTIONS  CALCULATED. 

C  RFR  =  "N"  =  REAL  PART  OF  INDEX  OF  REFRACTION. 

C  RFI  =  "K"  =  IMAGINARY  PART  OF  INDEX  OF  REFRACTION 
C  DISP  =  DISTANCE  OF  CENTER  OF  SPHERE  ABOVE  MIRROR. 

C  LMAX  =  NUMBER  OF  MODES  IN  CALCULATION. 

C  C2EXT  =  EXTINCTION  CROSS  SECTION  FOR  PARTICLE  ON  MIRROR. 

C  C2ABS  =  ABSORPTION  CROSS  SECTION  FOR  PARTICLE  ON  MIRROR. 

C  C2SCA  =  SCATTERING  CROSS  SECTION  FOR  PARTICLE  ON  MIRROR. 

C  CEXT  =  EXTINCTION  CROSS  SECTION  FOR  ISOLATED  PARTICLE. 

C  CABS  =  ABSORPTION  CROSS  SECTION  FOR  ISOLATED  PARTICLE. 

C  CSCA  =  SCATTERING  CROSS  SECTION  FOR  ISOLATED  PARTICLE. 

C  CDIF  =  DIFFERENTIAL  CROSS  SECTION  FOR  PARTICLE  ON  MIRROR. 

C  SDIF  =  DIFFERENTIAL  CROSS  SECTION  FOR  ISOLATED  PARTICLE. 

C  DEFAULT  VALUES  OF  INPUT  PARAMETERS 
WIr=l.  0 
RAD=1.0 
DELANG=5. 

NANG=19 

RFR=1.46 

RFI=0. 

DISP=RAD 

C  FIXED  INPUT  PARAMETERS 
999  LMAX=0 

ACRCY=1.E2 

LDIM=60 

C  READ  INPUT  VARIABLES 

CALL  INPVAR(WL, RAD, DELANG, NANG, RFR, RFI, DISP) 

IF (WL . LE . 0 . )  GO  TO  2000 
PI=ACOS (-1.0) 

KAY =2 . *PI/WL 
SZP=RAD*KAY 
AREA=PI*RAD*RAD 
CON=PI/180 . 

DO  50  1=1, NANG 

50  THETA(I)=FLOAT(I-l) *DELANG*CON 

IF (LMAX, LE . 0)  CALL  LMAXX(SZP,LDIM,ACRCY,LMAX) 

C  WRITE  THE  INPUT  VARIABLES 
WRITE (6,  65) 

65  FORMAT ( / 1  PROGRAM  MIRMX'/) 

WRITE (6, 100)  WL, RAD, SZP, DISP, LMAX, RFR, RFI 
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100  FORMAT ( 1  WL  =  ',1PE12.3/'  RAD  =  ’,1PE12.3/’  SZP  =  \1PE12.3/ 

+  •  DISP  =  ’ , 1PE11 . 3/ ’  LMAX  =  ',14/*  RFR  =  *,1PE12.3/ 

+  ’  RFI  =  ’ , 1PE12. 3/) 

CALL  ABMATR (LMAX, KAY, DISP) 

CALL  MIEUV 

CALL  GFMATX (LMAX, DISP, KAY, 1) 

CALL  FSOLVE (LMAX) 

CALL  CRSEC (SZP, DISP, KAY, LMAX) 

C  PRINT  THE  CROSS  SECTIONS 
Q2SCA=Q2EXT-Q2ABS 
C2EXT=AREA*Q2EXT 
C2SCA=AREA*Q2SCA 
C2ABS=AREA*Q2ABS 
WRITE (6, 600)  C2EXT,C2SCA,C2ABS 
600  FORMAT (//’  PARTICLE  ON  A  MIRROR  CROSS  SECTIONS’/ 

+’  C2EXT  =  1 , 1PE12 . 5, 5X, 1 C2SCA  =  ’ , 1PE12 . 5, 5X, • C2ABS  *  ’,1PE12.5/) 
CEXT=AREA*QEXT 
CSCA=AREA*QSCA 
CABS=AREA*QABS 

WRITE (6, 1050)  CEXT, CSCA, CABS 
1050  FORMAT (//’  SINGLE  PARTICLE  CROSS  SECTIONS’/ 

+’  CEXT  =  ’ , 1PE12 . 5, 5X, • CSCA  »  • , 1PE12 . 5, 5X, • CABS  =  ’,1PE12.5/) 

IF (NANG.EQ. 0)  GO  TO  800 
WRITE (6,  998) 

998  FORMAT (///34X, ’SINGLE’ /34X, ’PARTICLE’ ) 

WRITE (6, 1000) 

1000  FORMAT (5X, ’ANGLE’, 6X, ’CROSS  SEC ’, 9X, ’ CROSS  SEC’/) 

DO  1010  1=1, NANG 
CDIF=AREA*DIFF (I) 

SDIF=AREA*DIFS (I) 

THETD=THETA (I ) /CON 
WRITE (6, 1020)  THETD, CDIF, SDIF 
1020  FORMAT  (F10 . 3, 1PE15 . 3, 3X, 1PE15 . 3) 

1010  CONTINUE 
800  WRITE (6, 1030) 

1030  FORMAT (1H  /’  *****************/) 

GO  TO  999 
2000  STOP 
END 

SUBROUTINE  LMAXX (X, LDIM, ACRCY , LMAX) 

C  THIS  SUBROUTINE  CALCULATES  THE  MAXIMUM  VALUE  OF  L  TO  BE  USED  IN  MIE 
C  SCATTERING  FOR  A  PARTICLE  OF  SIZE  PARAMETER  X  AND  ACCURACY 
C  PARAMETER  ACRCY.  LMAX  CANNOT  EXCEED  LDIM. 

L=INT (X) 

ZR=(2*L+1)/X 
TZR=2./X 
PROD=l . 

RJL=1. 

100  L=L+1 

ZR=ZR+TZR 

RJL=ZR-1 . /RJL 

P  ROD  =PROD  *  RJL 

IF  (L .  GT .  LDIM)  GO  TO  200 

IF (PROD. LT. ACRCY)  GO  TO  100 

LMAX=L 

RETURN 

200  WRITE (6, 400) 

400  FORMAT (’  ERROR  -  LMAX  IS  LARGER  THAN  THE  DIMENSIONED  ARRAYS’/ 

+9X, ’ (CHANGE  LDIM  AND  THE  CORRESPONDING  ARRAY  DIMENSIONS) ’ / 

+9X, ’ (OR  -  CHECK  FOR  ERRORS  IN  THE  INPUT  VARIABLES) ’//) 

STOP 

END 

SUBROUTINE  ABMATR (LMAX, KAY, DISP) 

C  THIS  SUBROUTINE  CLACULATES  THE  TRANSFORMATION  MATRICES  A  AND  B  FOR 
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C  TRANSLATION  OF  THE  VECTOR  SPHERICAL  WAVE  FUNCTIONS  ALONG  THE  Z-AXIS. 
C  [SEE  APPENDIX  OF  IEEE  AP-19  PAGE  378  (1971) ) 

COMPLEX  A, B, RBH, KD, SUMA, SUMB, COEF, IM 
REAL  KAY 

DIMENSION  RBH (0:400) ,FZ (400) ,FM(4QQ) 

COMMON/AB/  A  (60,  60),  B(60,60) 

CQD®40N/PRM/XMIN,  XMAX,  NOJ 
IM=  (0. ,1.0) 

D-2 . *DISP 
LMX2=2*LMAX+1 
KD=KAY*ABS (D) 

CALL  RBESH (KD, LMX2, RBH) 

C  GENERATE  A  AND  B  MATRICES 
DO  100  LL=1, LMAX 
XLL=LL 

DO  100  LM=1,LMAX 
XIM=LM 

FLL=LL* (LL+1 ) 

FLM=LM* (LM+1) 

FFLM=FLL+FLM 
SUMA= ( 0 .  ,  0 . ) 

SUMB= (0. ,  0 . ) 

CALL  THREE J (XLL, XLM, 0 . , 0 . , FZ ) 

CALL  THREE J (XLL, XLM, l.,-l.,FM) 

LMIN1=INT (XMIN+0 . 001 ) -1 
DO  200  J=1 , NO J, 2 
L=LMIN1+J 
FL=L* (L+l) 

TL1=2*L+1 

COEF= (IM** (-L) ) *TL1 *RBH (L) 

SUMA=SUMA+COEF* (FFLM-FL) *FZ ( J) *FM  ( J) 

SUMB=SUMB+COEF*FZ ( J) *FM( J) 

200  CONTINUE 

CF=FLOAT (2*IM+1) /FLM*SQRT (FLL/FLM) /2 . 

COEFCF*  (IM**  (LM-LL)  ) 

A  (LL,  LM)  =-COEF*  SUMA 
B  (LL,  LM) =COEF*SUMB* (2 . *D) * (KAY*IM) 

100  CONTINUE 
RETURN 
END 

SUBROUTINE  MIEUV 

C  THIS  SUBROUTINE  EVALUATES  THE  MIE  COEFFICIENTS  U  AND  V 
C  (WHERE  U  AND  V  ARE  THE  COEFFICIENTS  THAT  ARE  USUALLY 
C  REFERED  TO  AS  THE  MIE  A  AND  B  COEFFICIENTS.  THEY 
C  ARE  DEFINED  AS  IN  STRATTON'S  "ELECTROMAGNETIC  THEORY". 

COMPLEX  U,V,A,B,Q,D,G,N,ONE,IMG,RJ,RH,Z 
COMMON/RMIE/RAD, RFR, RFI,  WL, KAY,  LMAX 
COMMON/OUE/U  (100)  ,V (100) 

DIMENSION  A(100) ,B(100) ,RJ(100) ,RH(100) 

REAL  KAY 
ONE=(1.,0. ) 

IMG=(0.,1.) 

C  CALCULATE  A  AND  B 

N=CMPLX (RFR,  RFI ) 

Z=N*KAY*RAD 

CALL  RATRBJ (Z,  LMAX,  RJ) 

DO  100  L=1 , LMAX 
D=(L+1)/Z-RJ(L) 

A(L)=D*N 
100  B(L)=D/N 

C  CALCULATE  U  AND  V  COEFFICIENTS 
Z=RAD*KAY 

CALL  RATRBJ (Z, LMAX,  RJ) 

CALL  RATRBH (Z, LMAX,  RH) 

Q=- (CSIN(Z) /Z-CCOS(Z) ) *CEXP(-IMG*Z) / (IMG/Z+ONE) 
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DO  500  L=l, LMAX 
D= (L+l) /Z-RJ(L) 

G=  (L+l)  /Z-RH  (L) 

U  (L)  =-Q*  (A  (L)  -D)  /  (A  (L)  -G) 

V (L)=~Q* (B (L) -D) / (B  (L) -G) 

500  Q=Q*RJ(L)/RH(L) 

RETURN 

END 

***************************************** 

SUBROUTINE  GFMATX (LMAX, D, KAY) 

C  CALCULATE  THE  SUPERMATRICES  G  AND  F  FROM  THE  SUBMATRICES  A, B  AND  U,V 
COMMON/AB/A { 60 , 60 ) ,  B < 60,  60 ) 

CCMMON/CMIE/U (100) ,V(100) 

CCMMON/GFMAT/G (120, 120)  ,F<120) 

COMPLEX  A,  B,  G,  F,  KIMD,  U,  V,  EPS,  EPSTAR,  EPE,  EME 

COMPLEX  UA,UB,VA,VB,RHO,IM 

REAL  KAY 

DO  100  1=1, LMAX 

I2=I+I24AX 

DO  100  J=1,LMAX 

J2=J+LMAX 

G(I,  J)=A(J,  I) 

G(I,  J2)=B(J,I) 

G(I2,  J)=B(J,I) 

100  G(I2,  J2)=A(J,I) 

SGN=1.0 

DO  200  1=1, LMAX 

SGN=-SGN 

I2=I+LMAX 

G (I, I)=G (I, I) +SGN/U (I) 

200  G(I2, I2)=G (12, 12) -SGN/V (I) 

C 

C  CALCULATE  F  VECTOR 
C 

IM=  (0. , 1 . ) 

KIMD=IM*KAY*D 
EPS=CEXP (-KIMD) 

EPSTAR=CONJG  (EPS) 

SGN=1.0 

DO  300  1=1, LMAX 
SGN=-SGN 

EME=EPS-SGN*EPSTAR 
EPE=EPS + SGN*EP  STAR 
FL=FLOAT (2*1+1) /FLOAT (I* (1+1) ) 

rHOFL*  (IM**(I+1) ) 

F  (I)  =EME*RHO 
300  F (I+LMAX) =EPE*RHO 
RETURN 
END 

SUBROUTINE  F SOLVE (LMAX) 

COMMON/GFMAT/G (120, 120)  ,F (120) 

COMPLEX  G,  F 
DIMENSION  WA(120) 

LMX2=2*LMAX 

CALL  LEQT3  C (G, LMX2, 120, F,  1, 120, 0, WA,  IER) 

RETURN 

END 

SUBROUTINE  CRSEC (X, DISP, KAY, LMAX) 

C  CALCULATE  CROSS  SECTIONS 

COMMON/GFMAT/G (120, 120) ,F(120) 

COMMON/D IFF /NANG, THETA { 400 ) , DIFF (400), DIFS (400),  QEXT,  QSCA,  QABS 
COMMON/CRS2/Q2EXT, Q2ABS 
CCX4MON/CMIE/U  (100)  ,V(100) 

COMPLEX  G,F, IM, IKD, ST1, ST2, SP1, SP2, ST,SP, Al,Bl, A2, B2,U, V, EXPS, EXMN 
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COMPLEX  IL, IFL, ZERO 
REAL  KAY 

CONST=l . 0/  (2 . *ACOS (-1 . ) *X*X) 

LMX2=2*tjMAX 

IM=(0.,1.) 

ZERO=(0. ,  0. ) 

IKD=IM*KAY*DISP 

C  CALCULATE  SINGLE  PARTICLE  TOTAL  CROSS  SECTIONS 
SUMEX=0 . 0 
SUMSC=0 . 0 
DO  300  L=1,LMAX 
FL=FLOAT (2*L+1) 

PEX=-FL* (REAL (U (L) ) +REAL (V (L) ) ) 

PSC=FL* (U (L) *CONJG (U (L) ) +V (L) *CONJG (V (L) ) ) 
SUMEX~SUMEX+PEX 
300  SUM  >SUMSC+PSC 
X2=2./(X*X) 

QEXT=X2  *  SUMEX 
QSCA=X2  *  SUMSC 
QABS=QEXT-QSCA 

C  CALCULATE  2-PARTICLE  EXTINCTION  CROSS  SECTION 
STl=ZERO 
ST2=ZERO 
SGN=1.0 

DO  400  L=1,LMAX 
SGN=-SGN 
A2=F (L) 

B2=F (L+LMAX) 

A1=-SGN*A2 

B1=SGN*B2 

IFL=FLOAT (L* (L+l) ) * ( (-IM) **L) 

ST1=ST1+IFL* (Al+Bl) 

400  ST2=ST2+IFL* (A2+B2) 

EXPS=CEXP (IKD ) 

EXMN=CONJG  <;EXPS) 

Q2EXT=X2  * AIMAG (EXPS*ST1+EXMN*ST2) 

C  CALCULATE  2-PARTICLE  ABSORPTION  CROSS  SECTION 
Q2ABS=0 . 

DO  500  L=l, LMAX 
L2=L+LMAX 

AS=F (L) *CONJG (F (L) ) 

BS=F (L2 ) * CONOG (F (L2 ) } 

US=U (L) *CONJG (U (L) ) 

VS=V (L) *CONJG (V (L) ) 

FIr=FLOAT  ( (L*  (L+l)  )  **2)  /FLOAT  (2*L+1) 

500  Q2ABS=Q2ABS+FL* (AS/US* (US+REAL  (U (L) ) ) 
++BS/VS* (VS+REAL (V(L) ) ) ) 

Q2ABS=-X2*Q2ABS 
IF (NANG.EQ. 0)  RETURN 
DO  100  1=1, NANG 
Z=COS (THETA ( I ) ) 

STl=ZERO 
SPl=ZERO 
ST2=ZERO 
SP2=ZERO 
ST=ZERO 
SP=ZERO 
PIM=0 . 0 
PI=1 . 0 
SGN=1 . 0 

DO  200  L=l, LMAX 
SGN=-SGN 
A2=F  (L) 

B2=F (L+LMAX) 

A1=-SGN*A2 

B1=SGN*B2 
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FL=L 

S=Z*PI 

T=S-PIM 

PIP=S+ (FL+1 . ) /FL*T 
TAU=FL*T-PIM 

FLT=FLOAT (2*L+1) /FLOAT (L* (L+l) ) 

IL= (-IM) **L 

ST1=ST1+IL* (A1*PI+B1*TAU) 

SP1=SP1+IL* (A1*TAU+B1*PI) 

ST2=ST2+IL* (A2*PI+B2*TAU) 

SP2=SP2+IL* (A2*TAU+B2*PI) 

ST=ST-FLT* (U (L) *PI+V (L) *TAU) 

SP=SP-FLT* (V (L) *PI+U (L) *TAU) 

PIM=PI 
200  PI=PIP 

ST1=IM*ST1 

ST2=IM*ST2 

SP1=-SP1 

SP2=-SP2 

DIFS (I)=CONST* (ST*CONJG (ST) +SP*CONJG (SP) ) 

EXPS=CEXP (IKD*Z) 

EXMN=CONJG (EXPS) 

ST=EXPS*ST1+EXMN*ST2 
SP=EXPS* SP1+EXMN* SP2 

100  DIFF  (I) CONST*  (ST*CONJG  (ST)  +SP*CONJG  (SP) ) 

RETURN 

END 

SUBROUTINE  RBESH (Z, LMAX, RBH) 

C  THIS  SUBROUTINE  COMPUTES  THE  SPHERICAL  HANKEL  FUNCTION  HI 
C  (SEPT, 25,1989) 

DIMENSION  RBH (0:1) 

COMPLEX  Z , RBH, RBH1 , RBH2 ,  U,  IM 

U-(1.,0.) 

IM= (0. , 1 . ) 

CALL  RATXBH(Z, LMAX, RBH) 

RBH1=- IM/ Z  *  CEXP (IM*Z) 

DO  100  L=0,LMAX 
RBH2=RBH1*RBH (L) 

RBH (L) =RBH1 
100  RBH1=RBH2 

RBH (LMAX+1 ) =RBH2 

RETURN 

END 

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 

SUBROUTINE  RATXBH (Z, LMAX,  RH) 

C  THIS  SUBROUTINE  COMPUTES  THE  RATIO (L+l) /L  OF  RICCATI 
C  HANKEL  FUNCTIONS  UP  TO  ORDER  LMAX  AND  COMPLEX  ARGUMENT 
C  Z  BY  UPWARD  RECURRANCE.  (SEPT.  22,  1989) 

COMPLEX  Z,  ZR, TZR  U, IM,RHL, RH 
DIMENSION  RH(0:1) 

U=(l. 0,0.0) 

IM=(0. 0,1.0) 

RHL=U/Z-IM 
RH  ( 0 ) =RHL 
ZR=U/Z 
TZR=2. *ZR 
L=0 

100  L=L+1 

ZR=ZR+TZR 
RHL=ZR-U / RHL 
RH(L)=RHL 

IF (L.LT.LMAX)  GO  TO  100 

RETURN 

END 
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SUBROUTINE  RATRBJ (Z, LMAX, RJ) 

C  THIS  SUBROUTINE  COMPUTES  THE  RATIO  (L+l)/L  OF  RICCATI 
C  BESSEL  FUNCTIONS  "J"  UP  TO  ORDER  LMAX  AND  COMPLES  ARGUMENT  Z 
C  BY  DOWNWARD  RECURRANCE.  (SEPT  22,1989) 

COMPLEX  Z,ZR,TZR,U,RJL,RJ, PROD,  ERR 
DIMENSION  RJ (1) 

U= (1.0, 0.0) 

C  CALCULATE  INITIAL  L  VALUE  FOR  DOWNWARD  RECURSION 
RJL=U 
ZR=U/Z 
TZR=2.*ZR 
PROD=U 

ZR= (2*LMAX+1) *ZR 
L=LMAX 
100  L=L+1 

ZR=ZR+TZR 
RJL=ZR-U / RJL 
PROD=PROD  *  R JL 

E=SQRT (REAL (PROD) **2+AIMAG (PROD) **2) 

IF  (L.GT. 5000)  GO  TO  200 
IF  (E.LT.1.E8)  GO  TO  100 
C  DOWNWARD  RECURSION 
RJL=U/ (ZR+TZR) 

300  RJL=ZR-U/RJL 
L=L-1 
ZR=ZR-TZR 

IF (L.LE.LMAX)  RJ(L)=U/RJL 
IF (L.GT. 1)  GO  TO  300 
RETURN 

200  WRITE (6, 400) 

400  FORMAT («  L  EXCEED  5000  IN  RATSBJ' ) 

STOP 

END 

SUBROUTINE  RATRBH (Z, LMAX, RH) 

C  THIS  SUBROUTINE  COMPUTES  THE  RATIO (L+l)/L  OF  RICCATI 
C  HANKEL  FUNCTIONS  UP  TO  ORDER  LMAX  AND  COMPLEX  ARGUMENT 
C  Z  BY  UPWARD  RECURRANCE.  (SEPT.  22,  1989) 

COMPLEX  Z,  ZR,  TZR  U, IM,RHL,RH 
DIMENSION  RH (1) 

U=(l. 0,0.0) 

IM=(0. 0,1.0) 

RHL=U/Z-IM 

ZR==U/Z 

TZR=2.*ZR 

L=0 

100  L=L+1 

ZR=ZR+TZR 
RHL=ZR-U/RHL 
RH (L) =RHL 

IF (L.LT.LMAX)  GO  TO  100 

RETURN 

END 

SUBROUTINE  THREEJ(XJ2,XJ3,XM2,XM3,F) 

THIS  SUBROUTINE  COMPUTES  THE  3-J  SYMBOLS.  IT  IS  BASED  ON 
THE  METHOD  OF  SCHULTEN  AND  GORDON  J.MATH  PHYS.  VOL. 16,  P.1961. 

DIMENSION  F (500) , S (500) ,T (500) 

CQMMON/PRM/XMIN, XMAX, NOJ 

TEST  FOR  BAD  INPUT  DATA 

TSTl=AMOD (XJ3+ABS (XM3) , 1 . ) 

TST2=AMOD (XJ2+ABS (XM2) ,  1 . ) 
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IF((TST1.GT.0.) .OR. (TST2.GT.0.))  GO  TO  990 
IF(ABS(XM3)  .GT.XJ3)  GO  TO  990 
IF (ABS (XM2 ) . GT . X J2 )  GO  TO  990 

EPS=l.E-25 
EPSLN=1 . E-10 
XM1=-  (XM2+XM3) 

XMAX=XJ2+XJ3 

XMIN=AMAX1  (ABS  (XJ3-XJ2)  ,ABS  (XMl) ) 

NOJ=INT (XMAX-XMIN+1 . +EPSLN) 

IF (NOJ.EQ. 1)  GO  TO  900 

Cl= (XJ2-XJ3) **2 

C2= (XJ2+XJ3+1 . ) **2 

C3= (XJ2* (XJ2+1 . ) -XJ3* (XJ3+1 . ) ) *XM1 

C4=XM1*XM1 

C5=XM3-XM2 

IF ( (XM2.EQ. 0 . ) .AND. (XM3.EQ. 0 . ) )  GO  TO  800 

IF( (XM2.EQ.XM3) .AND. (XJ2.EQ.XJ3) )  GO  TO  800 

JCZ=  (NOJ+D/2 

NL=JCZ-1 

NR=NOJ-JCZ 

UPWARD  ITERATION 

S(1)=0. 

NLP1=NL+1 
DO  100  1=1, NLP 1 
X J=XMIN+FLOAT ( I -1 ) 

XJS=XJ*XJ 

A=SQRT (ABS ( (XJS-C1) * (C2-XJS) * (XJS-C4) ) ) 

B=- (2 . *XJ+1 . ) * (C3-XJ* (XJ+1 . )  *C5) 

YJ=XJ+1 . 

YJS=YJ*YJ 

A1=SQRT (ABS ( (YJS-C1) * (C2- YJS) * (YJS-C4 ) ) ) 
DEN=B+Y J*A*S (I ) 

IF (DEN.EQ. 0. )  DEN=EPS 
IF (XJ.NE. 0. )  S (1+1) =-XJ*Al/DEN 
IF (XJ.EQ.  0 . )  S (1+1) =-Al/C5 
100  CONTINUE 

DOWNWARD  ITERATION 


T  (NOJ)  =0 . 

DO  200  1=1, NR 
XJ=XMAX-FLOAT (1-1) 

XJS=XJ*XJ 

K=NOJ+l-I 

A=SQRT (ABS { (XJS-C1) * (C2-XJS) * (XJS-C4) ) ) 
B=- (2 . *XJ+1 . ) * (C3-XJ* (XJ+1 . ) *C5) 

YJ=XJ+1. 

YJS=YJ*YJ 

A1=SQRT (ABS ( (YJS-C1) * (C2-YJS) * (YJS-C4) ) ) 
DEN=B+XJ*A1*T (K) 

IF (DEN.EQ.0. )  DEN=EPS 
200  T (K-l ) =-Y J*A/DEN 

CALCULATE  F 

F(JCZ)=1. 

IF (NL.EQ. 0)  GO  TO  350 
DO  300  1=1, NL 
J=JCZ-I+1 

300  F(J-1)=S(J)*F(J) 

350  DO  400  1=1, NR 
J=JCZ+I-1 

400  F ( J+l ) =T ( J) *F ( J) 
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NORMALIZE  THE  RESULTS 
SUM=0. 

DO  500  1=1, NO J 
XJ=XMIN+FLOAT (1-1 ) 

500  SUM=SUM+ (2 . *XJ+1 . ) *F (I) *F (I) 

SUM=1 . /SQRT (SUM) 

DO  600  1=1, NO J 
600  F  (I)  =SUM*F  (I) 

GO  TO  700 

SPECIAL  CASES:  (M1=M2=M3=0)  OR  <M2=M3  AND  J2=J3) 

800  DO  810  1=1,  NOJ 
810  F (I) =0 . 

F  (1)  =1 . 

SUM=2 .  *XMIN+1 . 

DO  820  1=3, NOJ, 2 
XJ=XMIN+FLOAT (1-1 ) 

XJM1=XJ-1 . 

XJS=XJ*XJ 
XJS1=XJM1*XJM1 

AJ=SQRT  (ABS  ( (XJS-C1) * (C2-XJS) * (XJS-C4 ) ) ) 

AJ1=SQRT (ABS ( (XJS1-C1) * (C2-XJS1) * (XJS1-C4) ) ) 

F (I) =-XJ*AJl/ (XJM1*AJ) *F (1-2) 

820  SUM  =  SUM+ (2. *XJ+1 . ) *F (I) *F (I) 

SUM=1 . /SQRT (SUM) 

DO  830  1=1, NOJ, 2 
830  F (I) =SUM*F (I) 

CO  TO  700 

SPECIAL  CASE  NOJ=l 

900  F  (1 )  =1 .  / SQRT  (2 .  *XMAX+1 . ) 

DETERMINE  THE  CORRECT  SIGN 

700  K=INT (XJ2-XJ3-XM1+EPSLN) 

Sl= (-1 . ) **K 
S2=SIGN(l.,F(NOJ) ) 

SG=S1*S2 
DO  710  1=1, NOJ 
IF (ABS (F (I) ) .LE.EPS)  F(I)=0. 

710  F (I ) =SG*F (I) 

RETURN 

990  WRITE (6,  995) 

995  FORMAT (//'  - BAD  INPUT  DATA - '/) 

STOP 
END 

SUBROUTINE  INPVAR (WL, RAD, DELANG, NANG, RFR, RFI, DISP) 

C  THIS  SUBROUTINE  READS  (OR  GENERATES)  THE  INPUT  VARIABLES. 
C 

C  READ  THE  INPUT  VARIABLES 
WRITE (6, 700) 

WRITE (6, 770) 

WRITE (6, 701)  WL 

READ (5,*)  WL 

IF (WL . LE . 0 . )  GO  TO  1000 

WRITE (6, 702)  RAD 

READ (5, *)  RAD 

DISP=RAD 

WRITE (6, 703) 

WRITE (6, 705)  DELANG 
READ (5,*)  DELANG 
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WRITE (6, 706)  NANG 
READ (5,*)  NANG 
WRITE (6/  707) 

WRITE (6, 708)  RFR 
READ (5, *)  RFR 
WRITE (6, 709)  RFI 
READ  (5, *)  RFI 
WRITE (6, 714) 

WRITE (6, 712)  DISP 
READ (5,*)  DISP 
WRITE (6, 710) 

770  FORMAT (4X, 'TYPE  A  COMMA  TO  REPEAT  THE  PREVIOUS  VALUE'/ 

+4X, 'IF  NO  RESPONSE, TYPE  A  SECOND  COMMA’) 

700  FORMAT (4X, 'ENTER  WAVELENGTH  AND  PARTICLE  RADIUS  IN  SAME  UNITS  (USE 
+  WL=0  TO  END) ’ ) 

701  FORMAT (F10 . 3, 8X, 'WL=' ) 

702  FORMAT (F10. 3, 8X, 'RAD=') 

703  FORMAT (4X, 'ENTER  INITIAL  SCATTERING  ANGLE (DEG) ,  ANGLE  INCREMENT (DE 
+G)  AND  NUMBER  OF ' /4X, 'ANGLES (MAX=1000) .  USE  NANG=0  FOR  NO  ANGULAR 
+ CALCULATIONS ' ) 

705  FORMAT (F10. 3, 8X, ’DELANG=') 

706  FORMAT (110, 8X, 'NANG=' ) 

707  FORMAT (4X, 'ENTER  INDEX  OF  REFRACTION (M=N+I*K) ' ) 

708  FORMAT (F10. 3, 8X, 'N=') 

709  FORMAT (F10. 3, BX, 'K=' ) 

712  FORMAT (F10 . 3, 8X, ' DISP= ' ) 

714  FORMAT (4X,  'ENTER  DISPLACEMENT  FROM  MIRROR') 

710  FORMAT (///) 

1000  RETURN 

END 
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PROGRAM  MIRMX 


WL  = 

1 . 000E+00 

RAD  = 

1. 000E+00 

SZP  = 

6. 283E+00 

DISP  = 

1 . 000E+00 

LMAX  = 

12 

RFR  = 

1. 460E+00 

RFI  = 

0 . 000E+00 

PARTICLE  ON  A  MIRROR  CROSS  SECTIONS 

C2EXT  =  1 . 79407E+01  C2SCA  =  1.79407E+01  C2ABS  =  -2.13531E-12 

SINGLE  PARTICLE  CROSS  SECTIONS 

CF.XT  =  9. 22655E+00  CSCA  =  9.22655E+00  CABS  =  -7.58960E-13 


ANGLE 

CROSS  SEC 

SINGLE 
PARTICLE 
CROSS  SEC 

.000 

8 . 686E+01 

2 . 484E+01 

5.000 

7 . 719E+01 

2 . 237E+01 

10.000 

5.333E+01 

1 . 615E+01 

15.000 

2. 731E+01 

8 . 984E+00 

20.000 

9.610E+00 

3.575E+00 

25.000 

2. 426E+00 

1 . 077E+00 

30.000 

9 . 241E-01 

8.918E-01 

35.000 

6.348E-01 

1 . 546E+00 

40.000 

7. 631E-01 

1.848E+00 

45.000 

9. 904E-01 

1 . 481E+00 

50.000 

5.704E-01 

8.184E-01 

55.000 

3. 518E-01 

3.527E-01 

60.000 

6. 999E-01 

2 . 627E-01 

65.000 

4 . 685E-01 

3. 819E-01 

70.000 

9. 596E-02 

4.515E-01 

75.000 

1 . 861E-01 

3.653E-01 

80.000 

6. 593E-02 

2 . 091E-01 

85.000 

2. 744E-01 

1 . 191E-01 

90.000 

6. 653E-01 

1.329E-01 
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TECHNOLOGY  OPERATIONS 


The  Aerospace  Corporation  functions  as  an  "architect-engineer"  for  national  security 
programs,  specializing  in  advanced  military  space  systems.  The  Corporation's  Technology 
Operations  supports  the  effective  and  timely  development  and  operation  of  national  security 
systems  through  scientific  research  and  the  application  of  advanced  technology.  Vital  to  the 
success  of  the  Corporation  is  the  technical  staffs  wide-ranging  expertise  and  its  ability  to  stay 
abreast  of  new  technological  developments  and  program  support  issues  associated  with  rapidly 
evolving  space  systems.  Contributing  capabilities  are  provided  by  these  individual  Technology 
Centers: 

Electronics  Technology  Center:  Microelectronics,  solid-state  device  physics,  VLSI 
reliability,  compound  semiconductors,  radiation  hardening,  data  storage 
technologies,  infrared  detector  devices  and  testing;  electro-optics,  quantum 
electronics,  solid-state  lasers,  optical  propagation  and  communications;  cw  and 
pulsed  chemical  laser  development,  optical  resonators,  beam  control,  atmospheric 
propagation,  and  laser  effects  and  countermeasures;  atomic  frequency  standards, 
applied  laser  spectroscopy,  laser  chemistry,  laser  optoelectronics,  phase  conjugation 
and  coherent  imaging,  solar  cell  physics,  battery  electrochemistry,  battery  testing  and 
evaluation. 

Mechanics  and  Materials  Technology  Center:  Evaluation  and  characterization  of 
new  materials:  metals,  alloys,  ceramics,  polymers  and  their  composites,  and  new 
forms  of  carbon;  development  and  analysis  of  thin  films  and  deposition  techniques; 
nondestructive  evaluation,  component  failure  analysis  and  reliability;  fracture 
mechanics  and  stress  corrosion;  development  and  evaluation  of  hardened 
components;  analysis  and  evaluation  of  materials  at  cryogenic  and  elevated 
temperatures;  launch  vehicle  and  reentry  fluid  mechanics,  heat  transfer  and  flight 
dynamics;  chemical  and  electric  propulsion;  spacecraft  structural  mechanics, 
spacecraft  survivability  and  vulnerability  assessment;  contamination,  thermal  and 
structural  control;  high  temperature  thermomechanics,  gas  kinetics  and  radiation; 
lubrication  and  surface  phenomena. 

Space  and  Environment  Technology  Center:  Magnetospheric,  auroral  and  cosmic 
ray  physics,  wave-particle  interactions,  magnetospheric  plasma  waves;  atmospheric 
and  ionospheric  physics,  density  and  composition  of  the  upper  atmosphere,  remote 
sensing  using  atmospheric  radiation;  solar  physics,  infrared  astronomy,  infrared 
signature  analysis;  effects  of  solar  activity,  magnetic  storms  and  nuclear  explosions 
on  the  earth's  atmosphere,  ionosphere  and  magnetosphere;  effects  of  electromagnetic 
and  particulate  radiations  on  space  systems;  space  instrumentation;  propellant 
chemistry,  chemical  dynamics,  environmental  chemistry,  trace  detection; 
atmospheric  chemical  reactions,  atmospheric  optics,  light  scattering,  state-specific 
chemical  reactions  and  radiative  signatures  of  missile  plumes,  and  sensor  out-of- 
field-of-view  rejection. 


